function x_next = prediction_model(x, u, P_load)
% =========================================================================
% prediction_model.m - (多尺度版)
%
% 版本说明:
% 1. 【核心修改】将模型分为两个时间尺度进行计算。
% 2. 在 dt 控制尺度上: 首先计算一个理论的目标功率，并据此一次性地
%    计算出电池的功率分配和下一个时刻的 SOC 状态。
% 3. 在 dt/N 物理仿真尺度上: 通过内部循环，精确模拟燃料电池的
%    温度、衰减和实际输出功率的动态演化过程。
% =========================================================================

    % 从外部工作区加载模型参数
    persistent BattParams FCParams SimParams
    if isempty(BattParams)
        BattParams = evalin('base', 'BattParams');
        FCParams   = evalin('base', 'FCParams');
        SimParams  = evalin('base', 'SimParams');
    end

    dt = SimParams.Ts;
    epsilon = 1e-9; % 定义一个小的正数用于数值稳定
    
    % --- 提取函数调用开始时的初始状态 ---
    SOC_k        = x(1);
    T_k          = x(2);
    alpha_decay_k = x(3);
    P_fc_k       = x(4);
    
    % --- 控制输入和扰动在整个dt周期内保持不变 ---
    q_H2_in = u(1);
    q_air   = u(2);
    q_W     = u(3);
    
    % =====================================================================
    % 步骤 1: 在 dt 尺度上，计算目标功率并更新 SOC
    % =====================================================================

    % 1a. 基于初始状态(T_k)计算理论上的目标净功率
    I_target = (2 * FCParams.u_H2_util * FCParams.F_faraday * q_H2_in) / FCParams.N_cells;
    q_v_target = q_air * (FCParams.R_gas * FCParams.T0) / FCParams.p_s;
    kappa = FCParams.kappa_air;
    p_ratio = FCParams.p_c / FCParams.p_s;
    P_ad_target = (kappa / (kappa - 1)) * FCParams.p_s * q_v_target * (p_ratio^((kappa - 1) / kappa) - 1);
    P_comp_target = P_ad_target / FCParams.eta_ad;
    
    T_c_target = T_k - 273.15;
    log_p_h2o_target = 2.95e-2*T_c_target - 9.18e-5*T_c_target^2 + 1.44e-7*T_c_target^3 - 2.18;
    p_H2O_target = 10^log_p_h2o_target;
    
    p_H2_denom_target = exp(1.635 * (I_target / FCParams.A_cell_area) / (T_k^1.334)) * (p_H2O_target / FCParams.p_a);
    p_O2_denom_target = exp(4.192 * (I_target / FCParams.A_cell_area) / (T_k^1.334)) * (p_H2O_target / FCParams.p_c);

    safe_margin = 1.0 - 1e-6;
    p_H2_denom_target = min(p_H2_denom_target, safe_margin);
    p_O2_denom_target = min(p_O2_denom_target, safe_margin);

    p_H2_target = 0.5 * p_H2O_target * ( (1 / p_H2_denom_target) - 1 );
    p_O2_target = p_H2O_target * ( (1 / p_O2_denom_target) - 1 );
    
    p_H2_target = max(p_H2_target, epsilon);
    p_O2_target = max(p_O2_target, epsilon);

    C_O2_target = p_O2_target / (5.08e6 * exp(-498 / T_k));
    C_O2_target = max(C_O2_target, epsilon);
    
    E_nernst_target = 1.229 - 0.85e-3 * (T_k - 298.15) + 4.3085e-5 * T_k * log(p_H2_target * p_O2_target^0.5 + epsilon);
    U_act_target = -(-0.9 + 0.0024*T_k + 7.3e-5*T_k*log(C_O2_target) - 0.00023*T_k*log(I_target + epsilon));
    U_ohmic_target = (0.01605 - 3.5e-5*T_k + 8e-5*I_target) * FCParams.l_membrane;
    J_target = I_target / FCParams.A_cell_area;
    ratio_J_target = min(J_target / FCParams.J_max, 0.9999);
    U_con_target = -FCParams.B_const * log(1 - ratio_J_target);
    
    U_cell_target = E_nernst_target - U_act_target - U_ohmic_target - U_con_target;
    U_stack_target = FCParams.N_cells * U_cell_target;
    P_fc_stack_target = U_stack_target * I_target;
    P_system_target = P_fc_stack_target - P_comp_target;
    
    % 1b. 根据目标功率，一次性计算出下一个时刻的 SOC
    P_bat = P_load - P_system_target;
    I_bat = P_bat / interp1(BattParams.OCV_SOC_map(1,:), BattParams.OCV_SOC_map(2,:), SOC_k, 'linear', 'extrap');
    soc_next = SOC_k - (I_bat * dt) / BattParams.Capacity_As;

    % =====================================================================
    % 步骤 2: 在 dt/N 尺度上，高分辨率模拟 FC 内部状态的演化
    % =====================================================================
    N=4;
    dt_small = dt / N;
    
    % 初始化内部状态
    T_current       = T_k;
    alpha_current   = alpha_decay_k;
    p_fc_current    = P_fc_k;
    
    for i = 1:N
        % 2a. 在每个子步长，都根据当前的内部温度 T_current 重新计算所有物理量
        I_sub = (2 * FCParams.u_H2_util * FCParams.F_faraday * q_H2_in) / FCParams.N_cells;
        P_comp_sub = P_comp_target; % 假设辅助功耗在dt内不变

        % ... (此处省略了与 1a 中完全相同的电化学模型计算，因为结果是一样的)
        % 为了代码简洁，直接复用 P_system_target 作为每个子步长的目标功率
        % 注意：一个更复杂的模型可能会让目标功率也随温度变化，但当前模型下它是固定的
        p_fc_target_sub = P_system_target;

        % 2b. 更新FC内部状态 (欧拉前向积分)
        
        % 状态4: 燃料电池实际净功率 (这里简化为瞬时达到子步长目标，因为没有tau)
        p_fc_new = p_fc_target_sub;
        
        % 状态2: 电堆温度 T_stack
        U_cell_sub = U_cell_target; % 复用目标计算时的单体电压
        delta_Q_stack_sub = FCParams.N_cells * I_sub * (1.481 - U_cell_sub);
        delta_Q_rad_sub = FCParams.delta_emissivity * FCParams.sigma_boltzmann * FCParams.a_rad_area * (T_current^4 - FCParams.T0^4);
        delta_Q_W_sub = q_W * FCParams.Cp_water * (T_current - FCParams.T0);
        delta_Q_net_sub = delta_Q_stack_sub - delta_Q_rad_sub - delta_Q_W_sub;
        T_new = T_current + (1/FCParams.C_thermal) * delta_Q_net_sub * dt_small;

        % 状态3: 燃料电池累积衰减 alpha_decay
        delta_p_sub = p_fc_new - p_fc_current;
        p_stack_for_decay_sub = P_fc_stack_target; % 复用目标计算时的堆栈功率
        delta_alpha_sub = (FCParams.C1 * p_stack_for_decay_sub^2 + FCParams.C2 * delta_p_sub^2) * dt_small;
        alpha_new = alpha_current + delta_alpha_sub;
        
        % 2c. 更新内部状态，为下一个子步长做准备
        T_current     = T_new;
        alpha_current = alpha_new;
        p_fc_current  = p_fc_new;
    end

    % =====================================================================
    % 步骤 3: 返回组合后的最终状态向量
    % =====================================================================
    x_next = [soc_next; T_current; alpha_current; p_fc_current];
end

